Three-stage ultrafast demagnetization dynamics in a monolayer ferromagnet

Intense laser pulses can be used to demagnetize a magnetic material on an extremely short timescale. While this ultrafast demagnetization offers the potential for new magneto-optical devices, it poses challenges in capturing coupled spin-electron and spin-lattice dynamics. In this article, we study the photoinduced ultrafast demagnetization of a prototype monolayer ferromagnet Fe3GeTe2 and resolve the three-stage demagnetization process characterized by an ultrafast and substantial demagnetization on a timescale of 100 fs, followed by light-induced coherent A1g phonon dynamics which is strongly coupled to the spin dynamics in the next 200–800 fs. In the third stage, chiral lattice vibrations driven by nonlinear phonon couplings, both in-plane and out-of-plane are produced, resulting in significant spin precession. Nonadiabatic effects are found to introduce considerable phonon hardening and suppress the spin-lattice couplings during demagnetization. Our results advance our understanding of dynamic charge-spin-lattice couplings in the ultrafast demagnetization and evidence angular momentum transfer between the phonon and spin degrees of freedom.


Section S.1 Time-dependent density functional theory
Section S. 4 The FFT analysis for the magnetization dynamics Figure S3.The FFT spectrum of the magnetization oscillation during 200-800 fs.

S.1. TIME-DEPENDENT DENSITY FUNCTIONAL THEORY
The time evolution of electron wave functions is governed by the time-dependent Kohn-Sham (TDKS) equation [1]: where velocity gauge is used and the external field appears in the kinetic term in the form of vector potential A().
The propagation of TDKS orbitals is implemented on the adiabatic basis  ,k (r, ) where  and n denote the TDKS band index and the basis index, respectively.k refers to the reciprocal momentum index and  ,k () the time dependent coefficients.The adiabatic basis is calculated on the fly at each ionic step by diagonalizing the Hamiltonian: where  ,k () is the eigenvalue.Thus, the charge density can be calculated with  ,k () and  ,k (r, ): where   represents the occupation number of time-dependent Kohn-Sham wavefunctions.
For ions that are much heavier than electrons, their motions are treated classically on an averaged potential energy surface determined by the electronic distribution according to the Ehrenfest theorem.The nuclear positions are updated following the Hellmann-Feynman theorem [2]: where   and R  are the mass and position of the th ion.Eq. (S1) and Eq.(S5) represent the coupled electron-ion motion.
The Gaussian-envelop laser pulse we applied has the following waveform.

S.2. TDDFT IN THE PLANE WAVE BASIS AND THE EVOLUTION OPERATOR
In the representation of plane wave (PW) basis sets {G}, the time-dependent Kohn-Sham states  ,k (G, ) at each k point can be expanded in the adiabatic basis sets { ,k (G,  1 )}, which are the eigenstates of Hamiltonian  k (G,  1 ) [3,4]: where the coefficient  ,k () =  ,k (G,  1 ) |  ,k (G, ) .Writing  ,k () in the form of the coefficient matrix  k (), the time-dependent Kohn-Sham (TDKS) equation in the adiabatic basis At an infinitesimal time interval (Δ =  2 −  1 ), the Hamiltonian in the adiabatic basis can be regarded as varying linearly with time Using the evolution operator  k ( 2 ,  1 ), the Eq.S7 can be written as According to Crank-Nicholson algorithm,  k (  + ,   ) can be expanded as with  ′ =   + /2.Under the condition of sufficiently small step ( ≪ 1), we interpolate   (in our case   = 1000) copies during the time range Δ =  2 −  1 , i.e.   =  1 +  and  = Δ/  .Now that  k ( 2 ,  1 ) can be expressed as: From Eq. S11, the time step of the ionic system is Δ while the time step of electronic systems is  = Δ/  .Therefore in our calculations, the timesteps for the electronic and ionic systems differ by three orders of magnitude.

S.3. THE SPIN-LATTICE COUPLING IN BULK FGT
To gain a better understanding of the strength of spin-lattice coupling in the layered structure, we perform static calculations to obtain the variation of the magnetic moment in the -direction (Δ  ) in monolayer and bulk FGT for the phonon modes of our interest.The  1g phonon modes in bulk FGT are demonstrated in Fig. S2, with the two layers having either in-phase or out-of-phase atomic motions.As shown in Fig. S2, the variation of magnetic moment upon the same phonon displacement of the in-phase  1g -1 phonon mode of bulk FGT is very similar to that observed in the monolayer (Fig. 3 in the main text).However, the spin-phonon coupling strength for the in-phase  1g -2 mode in bulk FGT is roughly one-fifth of that in the monolayer.This indicates that the interlayer interaction might change the electron screening in bulk metallic magnetic materials, which suppresses the spin-phonon coupling compared to that in the monolayer.This aligns with previous work where strong spin-phonon coupling is mostly observed in fewer-layer samples [5].
As for the out-of-phase  1g phonons, the variation of magnetization (Δ  ) on phonon displacement shows a quadratic dependence in bulk FGT, and the amplitude is at least one order of magnitude smaller, due to the cancellation effect between adjacent layers.Therefore, we believe the monolayer is an ideal platform to study the ultrafast spin dynamics and can reveal the fundamental physics involved, yet a full-scale TDDFT simulation on bulk and few-layer FGT is highly desirable that we leave for future investigations.

S.4. THE FFT ANALYSIS FOR THE MAGNETIZATION DYNAMICS
As Fig. S3 shows, the FFT spectrum of the magnetization dynamics after 200 fs has been examined, in which there are two dominant peaks with frequency approximately  peak1 = 4.72 THz and  peak2 = 8.89 THz.The two frequency peaks correspond to the occurrence of the two A 1g phonon modes, which highlights the dominant role of A 1g phonon modes in the dynamics of the second-stage demagnetization.Besides, the two peak frequencies are slightly increased compared with the phonon eigenfrequency in the ground state (4.58 THz and 8.87 THz).These results are consistent with the phonon hardening effect in Fig. 5 in the main text.

S.5. A COMPLETE PRESENTATION OF THE IN-PLANE PHONON DYNAMICS
To explicitly examine the in-plane lattice vibration, we show in Fig. S4 the time-dependent projected intensities for all phonon modes.More specifically, the in-plane Raman-active A 2g phonon modes (denoted by yellow and brown lines) are excited upon imposing the laser pulse.
During the after-pulse dynamics, two sets of doubly-degenerat infrared active phonons (in green and purple lines) are activated and their amplitudes increase over time, especially for the modes in purple after 800 fs.This is likely arising from phonon coupling with the initially activated Ramanactive modes with large amplitudes.Therefore, there exist the in-plane phonons throughout the process but the chiral ones which affect the magnetization dynamics start to play a dominant role in the demagnetization dynamics after 800 fs.

S.6. THE MAGNETIZATION DYNAMICS TO 1250 FS
Due to the substantial computational cost, we have carried out simulations of the temporal magnetization dynamics at a longer timescale to 1250 fs, shown in Fig. S5.We do observe a recovery of the -component of magnetic moment   after 1000 fs.This is mostly likely due to the thermalization of lattice vibrations after 1 ps, from which the chiral phonon, or the circular motion of the Fe atoms becomes incoherent.This is an indication that the chiral phonons driven from nonlinear phononics might not maintain a strong coherence and the phase is strongly coupled with other phonon modes.If one can selectively drive the chiral phonons as proposed in for example Ref. [6], the spin precession could be better manipulated.

S.7. ROLE OF ELECTRON EXCITATION IN THE SPIN-PHONON COUPLING
To understand the suppressed spin-phonon coupling in the excited states, we make a direct comparison by tracking the magnetization dynamics when electron excitation is absent.Specifically, we take the transient lattice structure for each time step and calculate the magnetization, but taking the electron occupation numbers as in their ground state.As can be seen in Fig. S6, the magnetization did not show an obvious demagnetization but rather a fluctuation near the initial   value.Therefore, we believe the high electron temperature, or more specifically, the nonequilibrium electron distribution will affect the spin-lattice coupling.

S.9. NONLINEAR PHONONICS UNDER DIFFERENT PUMP FLUENCES
In this section, we discuss in more detail the phonon and spin dynamics in the weak and intense excitation regimes.In the main text, we demonstrate distinctive magnetization dynamics under

Figure S1 .
Figure S1.Waveform of the applied electric field.Section S.2 TDDFT in the plane wave basis and the evolution operator Section S.3 The spin-lattice coupling in bulk FGT Figure S2.The spin-phonon interaction in bulk FGT.

Figure S4 .
Figure S4.The time-dependent in-plane phonon dynamics.Section S.6 The magnetization dynamics to 1250 fs Figure S5.The time-dependent magnetization dynamics for the time range 0-1250 fs.Section S.7 Role of electron excitation in the spin-phonon coupling Figure S6.The role of electron excitation.Section S.8 Magnetization dynamics contributed from different atoms Figure S7.The atom-resolved magnetization dynamics of FGT.Section S.9 Nonlinear phononics under different pump fluences Figure S8.The excited A 1g and E 2u phonon amplitudes with increasing laser fluence.

Figure S9 .
Figure S9.The different behaviors of in-plane phonons under weak and intense excitation.

FIG. S1 .
FIG. S1.Waveform of the applied electric field linearly polarized along the -axis of the crystal cell.

FIG. S3 .
FIG. S3.The FFT spectrum of the magnetization oscillation during 200-800 fs.Two dominant FFT peaks of  peak1 = 4.72 THz and  peak2 = 8.89 THz are corresponding to the two A 1g phonon modes.

FIG. S4 .
FIG. S4.The time-dependent in-plane phonon dynamics.The yellow and brown lines denote the Ramanactive in-plane A 2g phonon modes; the green and purple lines denote two sets of doubly-degenerate infrared in-plane E 2u phonon modes.

FIG. S5 .
FIG. S5.The time-dependent magnetization dynamics for the time range 0-1250 fs.The red, green, and blue lines represent   ,   , and   , respectively.

FIG. S6 .
FIG. S6.The role of electron excitation.Light-driven magnetization dynamics with electron excitation (in black lines) and w/o electron excitation (in red lines).

weak
FIG. S8.The excited A 1g phonon amplitudes with increasing laser fluence, with a quadratic dependence at the low intensity limit.The purple blocks show the amplitude of the infrared-active in-plane E 2u phonon modes excited from nonlinear phononics, showing a threshold behavior.The black dotted vertical line indicates the threshold for the pump fluence, while the purple dotted horizontal line schematically shows the threshold for the E 2u phonon amplitude to induce chiral phonons and spin precession to guide the eye.